Appendix A: A Tutorial of Python
In [ ]:
#If using Colab, you can only pull files from your Google Drive
#Run the following code to give Colab access
from google.colab import drive
drive.mount('/content/drive')
#File paths should then look like this:
file_path = '/content/drive/MyDrive/LinAlg/data/file.png'
#As opposed to your computer's directory, which might look like this:
file_path = '/Users/yourname/LinAlg/data/file.png'
In [ ]:
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import sympy as sp
import scipy.integrate as integrate
from numpy.linalg import inv, det, eig, solve, svd
import statsmodels.api as sm
from google.colab import drive
drive.mount('/content/drive')
Mounted at /content/drive
Basic Commands and Calculations
In [ ]:
result = 1 + 4
print(result)
result2 = 2 + np.pi/4 - 0.8
print(result2)
x = 1
y = 2
z = 4
t = 2 * x**y - z
print(t)
u = 2
v = 3
print(u + v)
#Sine function in radians
print(np.sin(u * v)) #u*v = 6, so np.sin(6)
#Temperature data in a list
tmax = [77, 72, 75, 73, 66, 64, 59]
print(tmax)
#Generate a sequence
seq1 = np.arange(1, 9)
seq2 = np.arange(8) + 1
seq3 = np.arange(1, 9, 1)
seq4 = np.linspace(1, 8, 8)
seq5 = np.linspace(1, 8, 8)
print(seq1, seq2, seq3, seq4, seq5)
#Define a function
def samfctn(x):
return x * x
print(samfctn(4))
def fctn2(x, y, z):
return x + y - z / 2
print(fctn2(1, 2, 3))
5 1.9853981633974482 -2 5 -0.27941549819892586 [77, 72, 75, 73, 66, 64, 59] [1 2 3 4 5 6 7 8] [1 2 3 4 5 6 7 8] [1 2 3 4 5 6 7 8] [1. 2. 3. 4. 5. 6. 7. 8.] [1. 2. 3. 4. 5. 6. 7. 8.] 16 1.5
Basic Plots
In [ ]:
#Plot temperature data
tmax = [77, 72, 75, 73, 66, 64, 59]
plt.plot(range(1, 8), tmax)
plt.xlabel('Day')
plt.ylabel('Temperature')
plt.title('Temperature Data')
plt.show()
#More plot examples
#Plot the curve of y = sin(x) from -pi to 2*pi
x_vals = np.linspace(-np.pi, 2 * np.pi, 400)
y_vals = np.sin(x_vals)
plt.plot(x_vals, y_vals)
plt.title('Plot of y = sin(x) from -pi to 2*pi')
plt.xlabel('x')
plt.ylabel('sin(x)')
plt.show()
#Define a function to plot
def square(x):
return x * x
x_vals2 = np.linspace(-3, 2, 400)
y_vals2 = square(x_vals2)
plt.plot(x_vals2, y_vals2)
plt.title('Plot of y = x^2')
plt.xlabel('x')
plt.ylabel('x^2')
plt.show()
#Plot a 3D surface
x = np.linspace(-1, 1, 100)
y = np.linspace(-1, 1, 100)
X, Y = np.meshgrid(x, y)
Z = 1 - X**2 - Y**2 # z = 1 - x^2 - y^2
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(X, Y, Z, cmap='viridis', edgecolor='none')
ax.view_init(elev=20, azim=330) # Set perspective angle
plt.title('3D Surface Plot')
plt.show()
#Contour plot
plt.contour(X, Y, Z)
plt.title('Contour Plot')
plt.xlabel('X')
plt.ylabel('Y')
plt.show()
#Filled contour plot (color map)
plt.contourf(X, Y, Z, cmap='viridis')
plt.title('Filled Contour Plot')
plt.xlabel('X')
plt.ylabel('Y')
plt.colorbar() #To add a color scale
plt.show()
Symbolic calculations
In [ ]:
#Derivative of x^2 with respect to x
x = sp.symbols('x')
expr = x**2
derivative = sp.diff(expr, x)
print(derivative) #Answer: 2x
#Assigning a function and differentiating it
fx = x**2
derivative_fx = sp.diff(fx, x)
print(derivative_fx) # Answer: 2x
#Change the expression and differentiate
fx2 = x**2 * sp.sin(x)
derivative_fx2 = sp.diff(fx2, x)
print(derivative_fx2) #Answer: 2x*sin(x) + x^2cos(x)
#Defining a function of two variables
y = sp.symbols('y')
fxy = x**2 + y**2
print(fxy) #Expression of the function in terms of x and y
#Partial derivatives
partial_x = sp.diff(fxy, x)
partial_y = sp.diff(fxy, y)
print(partial_x) #Answer: 2x
print(partial_y) #Answer: 2y
#Numerical integration of x^2 from 0 to 1
def square_func(x):
return x**2
integral_square, error_square = integrate.quad(square_func, 0, 1)
print(integral_square) #Answer: 1/3 or 0.3333333 with error details
#Numerical integration of cos(x) from 0 to pi/2
def cos_func(x):
return np.cos(x)
integral_cos, error_cos = integrate.quad(cos_func, 0, np.pi / 2)
print(integral_cos) #Answer: 1 with error details
2*x 2*x x**2*cos(x) + 2*x*sin(x) x**2 + y**2 2*x 2*y 0.33333333333333337 0.9999999999999999
Vectors and Matrices
In [ ]:
vec = np.array([1, 6, 3, np.pi, -3]) #4X1 column vector
#Sequences
seq1 = np.arange(2, 7) #sequence from 2 to 6
seq2 = np.arange(1, 11, 2) #sequence from 1 to 10 in incriments of 2
#Vector arithmetic
x = np.array([1, -1, 1, -1])
print(x + 1) #Add 1 to each element
print(2 * x) #Multiply each element by 2
print(x / 2) #Divide each element by 2
#Elementwise operations and dot product
y = np.arange(1, 5)
print(x * y) #Elementwise multiplication
dot_product = np.dot(x, y) #Dot product
print(dot_product)
#Transpose and matrix multiplication
xt = x.reshape(1, 4) #row vector
yt = y.reshape(4, 1) #column vector
print(np.dot(xt, y.reshape(-1,1))) #1x1 matrix (dot product)
print(np.dot(x.reshape(-1,1), y.reshape(1, -1))) #4x4 outer product
#Matrix creation
my = y.reshape(2, 2, order='F') #column-major order
print(my)
print(my.shape)
print(my.flatten(order='F')) #column-wise flattening
#Matrix operations
mx = np.array([[1, 1], [-1, -1]])
print(mx * my) #elementwise multiplication
print(mx / my) #elementwise division
print(mx - 2 * my) #subtraction
print(np.dot(mx, my)) #matrix multiplication
#Determinant
print(det(my))
#Inverse
myinv = inv(my)
print(myinv)
print(np.dot(myinv, my)) #should be identity
#Diagonal elements
print(np.diag(my))
#Eigen decomposition
eigvals, eigvecs = eig(my)
print("Eigenvalues:", eigvals)
print("Eigenvectors:\n", eigvecs)
#SVD
U, D, Vt = svd(my)
print("Singular values:", D)
print("Left singular vectors (U):\n", U)
print("Right singular vectors (V^T):\n", Vt)
#Solve linear system my * x = b
b = np.array([1, 3])
ysol = solve(my, b)
print("Solution:", ysol)
print("Verification:", np.dot(my, ysol)) #should match b
[2 0 2 0] [ 2 -2 2 -2] [ 0.5 -0.5 0.5 -0.5] [ 1 -2 3 -4] -2 [[-2]] [[ 1 2 3 4] [-1 -2 -3 -4] [ 1 2 3 4] [-1 -2 -3 -4]] [[1 3] [2 4]] (2, 2) [1 2 3 4] [[ 1 3] [-2 -4]] [[ 1. 0.33333333] [-0.5 -0.25 ]] [[-1 -5] [-5 -9]] [[ 3 7] [-3 -7]] -2.0 [[-2. 1.5] [ 1. -0.5]] [[1. 0.] [0. 1.]] [1 4] Eigenvalues: [-0.37228132 5.37228132] Eigenvectors: [[-0.90937671 -0.56576746] [ 0.41597356 -0.82456484]] Singular values: [5.4649857 0.36596619] Left singular vectors (U): [[-0.57604844 -0.81741556] [-0.81741556 0.57604844]] Right singular vectors (V^T): [[-0.40455358 -0.9145143 ] [ 0.9145143 -0.40455358]] Solution: [ 2.5 -0.5] Verification: [1. 3.]
Simple Statistics
In [ ]:
#Generate 10 normally distributed numbers
x = np.random.normal(size=10)
print("x =", x)
#Basic statistics
print("Mean:", np.mean(x))
print("Variance:", np.var(x, ddof=1)) #Sample variance
print("Standard deviation:", np.std(x, ddof=1)) #Sample sd
print("Median:", np.median(x))
print("Quantiles:", np.quantile(x, [0, 0.25, 0.5, 0.75, 1.0]))
print("Range (min, max):", (np.min(x), np.max(x)))
print("Max value:", np.max(x))
#Boxplot
plt.boxplot(x)
plt.title("Boxplot of x")
plt.show()
#Generate 1000 normally distributed numbers
w = np.random.normal(size=1000)
#Statistical summary of 12 normal random numbers
summary_data = np.random.normal(size=12)
print("Summary statistics:")
print("Min:", np.min(summary_data))
print("1st Quartile:", np.quantile(summary_data, 0.25))
print("Median:", np.median(summary_data))
print("Mean:", np.mean(summary_data))
print("3rd Quartile:", np.quantile(summary_data, 0.75))
print("Max:", np.max(summary_data))
#Histogram of 1000 normal random numbers
plt.hist(w, bins=30, edgecolor='black')
plt.title("Histogram of 1000 Normal Random Numbers")
plt.xlabel("Value")
plt.ylabel("Frequency")
plt.show()
x = [-1.0103165 0.28189177 -0.22670931 1.98921454 -0.80252838 -1.00267092 1.17126383 0.09510206 -0.20308394 -2.11110327] Mean: -0.18189401338845723 Variance: 1.3673359913456582 Standard deviation: 1.1693314292131458 Median: -0.21489662546233773 Quantiles: [-2.11110327 -0.95263528 -0.21489663 0.23519434 1.98921454] Range (min, max): (np.float64(-2.111103272951647), np.float64(1.989214537995136)) Max value: 1.989214537995136
Summary statistics: Min: -2.1175707171911644 1st Quartile: -0.38502958620425637 Median: 0.4107859658065656 Mean: 0.2959054735259114 3rd Quartile: 0.9192220862603779 Max: 1.8557167217986883
Linear regression and linear trend line
In [ ]:
#2007-2016 data of the global temperature anomalies
#Source: NOAAGlobalTemp data
t = np.arange(2007, 2017) # Years from 2007 to 2016
T = np.array([0.36, 0.30, 0.39, 0.46, 0.33, 0.38, 0.42, 0.50, 0.66, 0.70])
#Linear regression model: T ~ t
X = sm.add_constant(t) #Add intercept term
model = sm.OLS(T, X).fit()
print(model.summary())
#Get slope (temperature change rate per year)
slope = model.params[1]
print(f"Temperature change rate: {slope:.5f} deg C/year or {slope*10:.2f} deg C/decade")
#Plot
plt.plot(t, T, marker='o', linestyle='-', label="Observed Data")
plt.plot(t, model.predict(X), color='red', linewidth=2, label="Linear Trend")
plt.xlabel("Year")
plt.ylabel("Temperature [°C]")
plt.title("2007–2016 Global Temperature Anomalies\nand Their Linear Trend [0.37 °C/decade]")
plt.legend()
plt.grid(True)
plt.show()
OLS Regression Results
==============================================================================
Dep. Variable: y R-squared: 0.680
Model: OLS Adj. R-squared: 0.640
Method: Least Squares F-statistic: 17.02
Date: Sun, 11 May 2025 Prob (F-statistic): 0.00332
Time: 06:04:49 Log-Likelihood: 12.076
No. Observations: 10 AIC: -20.15
Df Residuals: 8 BIC: -19.55
Df Model: 1
Covariance Type: nonrobust
==============================================================================
coef std err t P>|t| [0.025 0.975]
------------------------------------------------------------------------------
const -73.4269 17.909 -4.100 0.003 -114.725 -32.129
x1 0.0367 0.009 4.125 0.003 0.016 0.057
==============================================================================
Omnibus: 4.518 Durbin-Watson: 1.116
Prob(Omnibus): 0.104 Jarque-Bera (JB): 1.168
Skew: -0.147 Prob(JB): 0.558
Kurtosis: 1.352 Cond. No. 1.41e+06
==============================================================================
Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 1.41e+06. This might indicate that there are
strong multicollinearity or other numerical problems.
Temperature change rate: 0.03673 deg C/year or 0.37 deg C/decade
/usr/local/lib/python3.11/dist-packages/scipy/stats/_axis_nan_policy.py:430: UserWarning: `kurtosistest` p-value may be inaccurate with fewer than 20 observations; only n=10 observations were given. return hypotest_fun_in(*args, **kwds)
